Introduction to Bayesian Disease Measurement for Health Scientists - Day 1 | 2

CA18208 HARMONY Larissa Training School 2024 - https://harmony-net.eu/

Eleftherios Meletis
Polychronis Kostoulas
Evangelia Anyfanti

2024-04-22 & 2024-04-23

Bayesian Statistics

In this session we’ll see how we can estimate a probability of interest but in a Bayesian framework, i.e. using Bayes theorem.

Bayes’ theorem - P(A|B) = P(B|A)*P(A)/P(B)

Components

What we usually see/use

Where:

Bayesian Inference - Summary & Example

To estimate the posterior distribution P(theta|y) we need to:

Example 1

Hemophilia is a disease that exhibits X-chromosome-linked recessive inheritance. Consider a woman who has an affected brother, which implies that her mother must be a carrier of the hemophilia gene. We are also told that her father is not affected. Suppose she has two sons, neither of whom is affected.

Exercise 1: Bayesian apparent prevalence (ap) estimation

y out of n individuals test positive. Estimate the apparent prevalence in a Bayesian framework.

Let’s write our first JAGS model

ap_model <- 
  'model {
  
  # Define likelihood distribution of the data
  # JAGS Binomial distribution Arguments: p, n 
  
  y1 ~ dbin(ap1,n1)
  y2 ~ dbin(ap2,n2)
  
  # Specify prior distribution for par of interest 
  # Uniform (non-informative) prior distribution 
  ap1 ~ dbeta(1, 100000)
  ap2 ~ dbeta(1, 100000)


  #data# n1, y1, n2, y2
  #monitor# ap1, ap2
  #inits# ap1
  }
  '

Let’s run our first JAGS model

# Call JAGS
library(runjags)

# Provide Data 
n1 = 4000
y1 = 1200

n2 = 40
y2 = 12

# Initial values for par of interest

ap1 <- list(chain1=0.05, chain2=0.95)


# Run the model
results <- run.jags(ap_model, n.chains=1, 
                    burnin=5000, sample=10000)
## Warning: Convergence cannot be assessed with only 1 chain
# View results

# Plot results
plot(results)

# Print results
summary(results)
##          Lower95       Median      Upper95         Mean           SD Mode
## ap1 1.091972e-02 0.0115451934 0.0122377064 0.0115501328 3.355818e-04   NA
## ap2 6.553641e-05 0.0001266365 0.0002000996 0.0001296519 3.512836e-05   NA
##            MCerr MC%ofSD SSeff        AC.10 psrf
## ap1 4.270690e-06     1.3  6174  0.001341437   NA
## ap2 4.682192e-07     1.3  5629 -0.002783184   NA
# Convert the results to 'mcmc.list' object
mcmc_ap <- as.mcmc.list(results)

# Convert 'mcmc.list' to data frame for ggplot
df_ap <- as.data.frame(as.mcmc(mcmc_ap))

library(ggplot2)
library(tidyr)  # For pivot_longer

# Assuming df_ap contains the columns for ap1 and ap2
# If your data frame isn't set up this way, adjust the data extraction part accordingly.

# Convert the dataframe to long format
df_long <- pivot_longer(df_ap, cols = c(ap1, ap2), names_to = "parameter", values_to = "value")

library(ggplot2)

# Assuming df_long is already created and contains the parameters 'ap1' and 'ap2' in long format

# Generate a sample from a beta distribution
beta_data <- data.frame(value = rbeta(10000, 1, 10000))
beta_data$parameter <- "Beta Distribution"

# Combine this with your existing data
df_long <- rbind(df_long, beta_data)

# Plotting
p <- ggplot(df_long, aes(x = value, fill = parameter)) +
  geom_density(alpha = 0.6, adjust = 1) +  # Adjust can be tweaked for smoother curves
  scale_fill_manual(values = c(ap1 = "blue", ap2 = "red", `Beta Distribution` = "green")) +
  labs(title = "Overlay of Density Plots for Apparent Prevalences ap1, ap2, and Beta Distribution",
       x = "Prevalence Estimate",
       y = "Density") +
  theme_minimal()

# Print the plot
print(p)

Exercise 1 - Practical

Exercise 2: Bayesian true prevalence (tp) estimation

Assuming the absence of a perfect test we do not know how many individuals are truly positive/negative.

Instead we know that n individuals are tested with an imperfect test and y have a positive result.

Create a JAGS model for true prevalence estimation

Prior distributions

Data: n tested, y positive

Write JAGS model

tp_model <- 
  'model {
  y ~ dbin(ap,n)
  ap <- tp*Se + (1-tp)*(1-Sp)
  
  # Uniform (non-informative) prior distribution 
  tp ~ dbeta(1,1)
  # Informative priors for Se and Sp
  Se ~ dbeta(25.4, 3.4)
  Sp ~ dbeta(95, 5)
  
  #data# n, y
  #monitor# tp, Se, Sp
  #inits# tp, Se, Sp
  }
  '

Let’s run our JAGS model

# Call JAGS
library(runjags)
# Provide Data 
n = 4072
y = 1210
# Initial values for pars of interest
tp <- list(chain1=0.05, chain2=0.95)
Se <- list(chain1=0.05, chain2=0.95)
Sp <- list(chain1=0.05, chain2=0.95)

# Run the model
results <- run.jags(tp_model, n.chains=2, 
                    burnin=5000, sample=10000)
## View results
# Plot results
plot(results) # Click backwards to view all plots

# Print results
summary(results)
##      Lower95    Median   Upper95      Mean         SD Mode        MCerr MC%ofSD
## tp 0.2403075 0.2978380 0.3684788 0.3001641 0.03217121   NA 0.0012618569     3.9
## Se 0.7566369 0.8866616 0.9812408 0.8770171 0.06178538   NA 0.0022247881     3.6
## Sp 0.9076471 0.9524757 0.9866954 0.9497712 0.02142762   NA 0.0006325174     3.0
##    SSeff     AC.10     psrf
## tp   650 0.5536472 1.014044
## Se   771 0.4485606 1.002960
## Sp  1148 0.3088396 1.003026

Exercise 2 - Practical

Ok for today? - Any questions?

Sensitivity - Specificity estimation without a gold standard

Hui-Walter explained

https://www.youtube.com/watch?v=z6devQmW2xE&ab_channel=PolychronisKostoulas#

We will use the data/observations from the manuscript published back in 1980.

Hui-Walter (1980) dataset Table 1

pop_1 = matrix(nrow=3,ncol=3)
rownames(pop_1) = c("Mantoux_Test_Pos", "Mantoux_Test_Neg", "Total")
colnames(pop_1) = c("Tine_Test_Pos", "Tine_Test_Neg", "Total")

pop_1[1,1] = 14
pop_1[1,2] = 4
pop_1[2,1] = 9
pop_1[2,2] = 528
#Total rows and columns
pop_1[1,3] = pop_1[1,1] + pop_1[1,2]
pop_1[2,3] = pop_1[2,1] + pop_1[2,2]
pop_1[3,1] = pop_1[1,1] + pop_1[2,1]
pop_1[3,2] = pop_1[1,2] + pop_1[2,2]
N_1 = sum(pop_1[1,1] + pop_1[1,2] + pop_1[2,1] + pop_1[2,2])
pop_1[3,3] = N_1
pop_1
##                  Tine_Test_Pos Tine_Test_Neg Total
## Mantoux_Test_Pos            14             4    18
## Mantoux_Test_Neg             9           528   537
## Total                       23           532   555
## Now let's do pop_2
pop_2 = matrix(nrow=3,ncol=3)
rownames(pop_2) = c("Mantoux_Test_Pos", "Mantoux_Test_Neg", "Total")
colnames(pop_2) = c("Tine_Test_Pos", "Tine_Test_Neg", "Total")

pop_2[1,1] = 887
pop_2[1,2] = 31
pop_2[2,1] = 37
pop_2[2,2] = 367
#Total rows and columns
pop_2[1,3] = pop_2[1,1] + pop_2[1,2]
pop_2[2,3] = pop_2[2,1] + pop_2[2,2]
pop_2[3,1] = pop_2[1,1] + pop_2[2,1]
pop_2[3,2] = pop_2[1,2] + pop_2[2,2]
N_2 = sum(pop_2[1,1] + pop_2[1,2] + pop_2[2,1] + pop_2[2,2])
pop_2[3,3] = N_2
pop_2
##                  Tine_Test_Pos Tine_Test_Neg Total
## Mantoux_Test_Pos           887            31   918
## Mantoux_Test_Neg            37           367   404
## Total                      924           398  1322

Hui-Walter model

Model Specification (‘hw_definition’)

hw_definition <- c("model{
  Population_1 ~ dmulti(prob_1, N_1)
  Population_2 ~ dmulti(prob_2, N_2)
  
  #Population_1
  
  # Test1+ Test2+
    prob_1[1] <- (prev[1] * ((se[1])*(se[2]))) + ((1-prev[1]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_1[2] <- (prev[1] * ((se[1])*(1-se[2]))) + ((1-prev[1]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_1[3] <- (prev[1] * ((1-se[1])*(se[2]))) + ((1-prev[1]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_1[4] <- (prev[1] * ((1-se[1])*(1-se[2]))) + ((1-prev[1]) * ((sp[1])*(sp[2])))
    
    #Population_2
  
  # Test1+ Test2+
    prob_2[1] <- (prev[2] * ((se[1])*(se[2]))) + ((1-prev[2]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_2[2] <- (prev[2] * ((se[1])*(1-se[2]))) + ((1-prev[2]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_2[3] <- (prev[2] * ((1-se[1])*(se[2]))) + ((1-prev[2]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_2[4] <- (prev[2] * ((1-se[1])*(1-se[2]))) + ((1-prev[2]) * ((sp[1])*(sp[2])))

  prev[1] ~ dbeta(1, 1)
  prev[2] ~ dbeta(1, 1)
  
  se[1] ~ dbeta(1, 1)I(1-sp[1], )
  sp[1] ~ dbeta(1, 1)
  se[2] ~ dbeta(1, 1)I(1-sp[2], )
  sp[2] ~ dbeta(1, 1)

  #data# Population_1, Population_2, N_1, N_2
  #monitor# prev, prob_1, prob_2, se, sp
  #inits# prev, se, sp
  }
  ")

Let’s run the model

library('runjags')

Population_1 <- as.numeric(pop_1[1:2,1:2])
Population_2 <- as.numeric(pop_2[1:2,1:2])


prev <- list(chain1=c(0.05,0.99), chain2=c(0.95,0.05))
se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

results <- run.jags(hw_definition, n.chains=2)

Remember to check convergence and effective sample size!

pt <- plot(results)

pt$`prev[1].plot1`

pt$`prev[1].plot3`

print(pt[["prev[1].plot1"]])

print(pt[["se[1].plot1"]])

print(pt[["sp[1].plot1"]])

print(pt[["sp[1].plot3"]])

summary(results)
##               Lower95      Median    Upper95        Mean          SD Mode
## prev[1]   0.015134996 0.028044913 0.04307840 0.028627279 0.007283481   NA
## prev[2]   0.690915330 0.715582840 0.74133190 0.715430590 0.012890610   NA
## prob_1[1] 0.014157200 0.026360418 0.04027954 0.026920260 0.006791926   NA
## prob_1[2] 0.007735796 0.017439528 0.02914599 0.018038765 0.005652170   NA
## prob_1[3] 0.002429108 0.008405407 0.01684339 0.008968083 0.003920497   NA
## prob_1[4] 0.927018653 0.946654656 0.96376544 0.946072892 0.009512237   NA
## prob_2[1] 0.643579183 0.669379993 0.69426946 0.669237567 0.012986677   NA
## prob_2[2] 0.019577493 0.028400701 0.03778528 0.028660257 0.004672813   NA
## prob_2[3] 0.015911645 0.023935304 0.03248369 0.024196034 0.004240767   NA
## prob_2[4] 0.253458655 0.277790209 0.30180763 0.277906142 0.012333790   NA
## se[1]     0.955584310 0.968754974 0.98063872 0.968429806 0.006411650   NA
## se[2]     0.951673572 0.966191131 0.98003042 0.965888714 0.007232896   NA
## sp[1]     0.970809468 0.982881713 0.99337004 0.982256434 0.005949219   NA
## sp[2]     0.983163647 0.992115396 0.99852033 0.991520346 0.004173473   NA
##                  MCerr MC%ofSD SSeff         AC.10      psrf
## prev[1]   6.671775e-05     0.9 11918 -0.0141690806 1.0001733
## prev[2]   1.216458e-04     0.9 11229 -0.0055674909 1.0000695
## prob_1[1] 6.193725e-05     0.9 12025 -0.0148276346 1.0001792
## prob_1[2] 6.160316e-05     1.1  8418 -0.0028020366 0.9999652
## prob_1[3] 4.231125e-05     1.1  8586  0.0024832450 1.0000058
## prob_1[4] 9.451507e-05     1.0 10129 -0.0185115477 1.0000470
## prob_2[1] 1.144714e-04     0.9 12871  0.0008074865 1.0000947
## prob_2[2] 4.372839e-05     0.9 11419 -0.0020506095 1.0000674
## prob_2[3] 3.980346e-05     0.9 11351  0.0143988167 1.0000089
## prob_2[4] 1.124378e-04     0.9 12033 -0.0017992025 1.0001151
## se[1]     6.565136e-05     1.0  9538  0.0190230422 0.9999648
## se[2]     7.659965e-05     1.1  8916 -0.0020043666 1.0000175
## sp[1]     6.558031e-05     1.1  8229 -0.0026330896 0.9999660
## sp[2]     4.569984e-05     1.1  8340  0.0043092761 0.9999978

Exercise 1

Solution - Model Specification (‘hw_definition’)

hw_definition_1 <- c("model{
  Population_1 ~ dmulti(prob_1, N_1)
  Population_2 ~ dmulti(prob_2, N_2)
  
  #Population_1
  
  # Test1+ Test2+
    prob_1[1] <- (prev[1] * ((se[1])*(se[2]))) + ((1-prev[1]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_1[2] <- (prev[1] * ((se[1])*(1-se[2]))) + ((1-prev[1]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_1[3] <- (prev[1] * ((1-se[1])*(se[2]))) + ((1-prev[1]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_1[4] <- (prev[1] * ((1-se[1])*(1-se[2]))) + ((1-prev[1]) * ((sp[1])*(sp[2])))
    
    #Population_2
  
  # Test1+ Test2+
    prob_2[1] <- (prev[2] * ((se[1])*(se[2]))) + ((1-prev[2]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_2[2] <- (prev[2] * ((se[1])*(1-se[2]))) + ((1-prev[2]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_2[3] <- (prev[2] * ((1-se[1])*(se[2]))) + ((1-prev[2]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_2[4] <- (prev[2] * ((1-se[1])*(1-se[2]))) + ((1-prev[2]) * ((sp[1])*(sp[2])))

  prev[1] ~ dbeta(1, 1)
  prev[2] ~ dbeta(1, 1)
  
  se[1] ~ dbeta(5, 1)I(1-sp[1], )
  sp[1] ~ dbeta(5, 1)
  se[2] ~ dbeta(1, 1)I(1-sp[2], )
  sp[2] ~ dbeta(1, 1)

  #data# Population_1, Population_2, N_1, N_2
  #monitor# prev, prob_1, prob_2, se, sp
  #inits# prev, se, sp
  }
  ")

Population_1 <- as.numeric(pop_1[1:2,1:2])
Population_2 <- as.numeric(pop_2[1:2,1:2])


prev <- list(chain1=c(0.05,0.99), chain2=c(0.95,0.05))
se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

results_1 <- run.jags(hw_definition_1, n.chains=2)
# Remember to check convergence and effective sample size!
plot(results_1)  # Click backwards to view all plots

pt_1 <- plot(results_1)

pt_1$`prev[1].plot1`

pt_1$`prev[1].plot3`

print(pt_1[["prev[1].plot1"]])

print(pt_1[["se[1].plot1"]])

print(pt_1[["sp[1].plot1"]])

print(pt_1[["sp[1].plot3"]])

summary(results_1)
##               Lower95      Median    Upper95        Mean          SD Mode
## prev[1]   0.014826923 0.028134541 0.04306610 0.028749402 0.007438955   NA
## prev[2]   0.689752755 0.715504409 0.74021997 0.715359756 0.012947617   NA
## prob_1[1] 0.014773174 0.026470000 0.04116364 0.027035759 0.006934728   NA
## prob_1[2] 0.007476164 0.017319275 0.02904750 0.017931072 0.005663263   NA
## prob_1[3] 0.002048212 0.008427675 0.01675232 0.009013356 0.004000690   NA
## prob_1[4] 0.927204204 0.946435070 0.96494019 0.946019814 0.009646584   NA
## prob_2[1] 0.643645747 0.669390102 0.69435104 0.669222470 0.012932065   NA
## prob_2[2] 0.020076707 0.028527605 0.03795316 0.028734909 0.004622913   NA
## prob_2[3] 0.016287730 0.023785364 0.03256279 0.024051480 0.004207159   NA
## prob_2[4] 0.253883200 0.277833017 0.30226045 0.277991141 0.012376384   NA
## se[1]     0.955465800 0.969014218 0.98037901 0.968654657 0.006363729   NA
## se[2]     0.951839701 0.965934467 0.97960843 0.965740429 0.007153362   NA
## sp[1]     0.970736268 0.983012684 0.99345855 0.982375180 0.005961492   NA
## sp[2]     0.983262065 0.992083439 0.99895634 0.991469690 0.004255582   NA
##                  MCerr MC%ofSD SSeff         AC.10      psrf
## prev[1]   7.102587e-05     1.0 10970  0.0010270420 1.0000487
## prev[2]   1.220861e-04     0.9 11247  0.0081960861 1.0000494
## prob_1[1] 6.592711e-05     1.0 11064  0.0003184868 1.0000450
## prob_1[2] 6.215152e-05     1.1  8303  0.0154585463 0.9999592
## prob_1[3] 4.446067e-05     1.1  8097 -0.0029011841 1.0000610
## prob_1[4] 9.457490e-05     1.0 10404 -0.0097024767 1.0002971
## prob_2[1] 1.124436e-04     0.9 13227 -0.0012019476 1.0004216
## prob_2[2] 4.219356e-05     0.9 12004  0.0099645710 1.0003116
## prob_2[3] 4.173170e-05     1.0 10164 -0.0057853883 1.0004645
## prob_2[4] 1.117283e-04     0.9 12270  0.0046014984 1.0000558
## se[1]     6.870037e-05     1.1  8580 -0.0066740779 1.0005514
## se[2]     7.433448e-05     1.0  9261  0.0094238536 1.0002994
## sp[1]     6.615695e-05     1.1  8120  0.0153886562 0.9999518
## sp[2]     4.798811e-05     1.1  7864 -0.0019992100 1.0001384

Exercise 2

Solution Model Specification (‘hw_definition’)

hw_definition_2 <- c("model{
  Population_1 ~ dmulti(prob_1, N_1)
  Population_2 ~ dmulti(prob_2, N_2)
  
  #Population_1
  
  # Test1+ Test2+
    prob_1[1] <- (prev[1] * ((se[1])*(se[2]))) + ((1-prev[1]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_1[2] <- (prev[1] * ((se[1])*(1-se[2]))) + ((1-prev[1]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_1[3] <- (prev[1] * ((1-se[1])*(se[2]))) + ((1-prev[1]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_1[4] <- (prev[1] * ((1-se[1])*(1-se[2]))) + ((1-prev[1]) * ((sp[1])*(sp[2])))
    
    #Population_2
  
  # Test1+ Test2+
    prob_2[1] <- (prev[2] * ((se[1])*(se[2]))) + ((1-prev[2]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_2[2] <- (prev[2] * ((se[1])*(1-se[2]))) + ((1-prev[2]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_2[3] <- (prev[2] * ((1-se[1])*(se[2]))) + ((1-prev[2]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_2[4] <- (prev[2] * ((1-se[1])*(1-se[2]))) + ((1-prev[2]) * ((sp[1])*(sp[2])))

  prev[1] ~ dbeta(1, 1)
  prev[2] ~ dbeta(1, 1)
  
  se[1] ~ dbeta(1, 1)
  sp[1] ~ dbeta(1, 1)
  se[2] ~ dbeta(1, 1)
  sp[2] ~ dbeta(1, 1)

  #data# Population_1, Population_2, N_1, N_2
  #monitor# prev, prob_1, prob_2, se, sp
  #inits# prev, se, sp
  }
  ")

Population_1 <- as.numeric(pop_1[1:2,1:2])
Population_2 <- as.numeric(pop_2[1:2,1:2])


prev <- list(chain1=c(0.05,0.99), chain2=c(0.95,0.05))
se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

results_2 <- run.jags(hw_definition_2, n.chains=2)
# Remember to check convergence and effective sample size!
plot(results_2)  # Click backwards to view all plots

pt_2 <- plot(results_2)

pt_2$`prev[1].plot1`

pt_2$`prev[1].plot3`

print(pt_2[["prev[1].plot1"]])

print(pt_2[["se[1].plot1"]])

print(pt_2[["sp[1].plot1"]])

print(pt_2[["sp[1].plot3"]])

print(pt_2[["se[2].plot1"]])

summary(results_2)
##               Lower95     Median    Upper95        Mean          SD Mode
## prev[1]   0.015011314 0.02810646 0.04296564 0.028699435 0.007325016   NA
## prev[2]   0.690992901 0.71563005 0.74209896 0.715575887 0.013008858   NA
## prob_1[1] 0.014169211 0.02643493 0.04025525 0.026983627 0.006826925   NA
## prob_1[2] 0.007875020 0.01739199 0.02930559 0.017992973 0.005620845   NA
## prob_1[3] 0.002184289 0.00832246 0.01673989 0.008913238 0.003969325   NA
## prob_1[4] 0.927010407 0.94661370 0.96400615 0.946110162 0.009511643   NA
## prob_2[1] 0.643445447 0.66943086 0.69458122 0.669312627 0.013060432   NA
## prob_2[2] 0.019986687 0.02843649 0.03768666 0.028651309 0.004561164   NA
## prob_2[3] 0.016203085 0.02396447 0.03282358 0.024238713 0.004293344   NA
## prob_2[4] 0.253234837 0.27767702 0.30196209 0.277797351 0.012439240   NA
## se[1]     0.955313036 0.96868036 0.98048564 0.968349287 0.006465680   NA
## se[2]     0.952246961 0.96611973 0.97977093 0.965881641 0.007055120   NA
## sp[1]     0.970243662 0.98293571 0.99280968 0.982307275 0.005915084   NA
## sp[2]     0.983237189 0.99219836 0.99874894 0.991581913 0.004223711   NA
##                  MCerr MC%ofSD SSeff         AC.10      psrf
## prev[1]   7.051710e-05     1.0 10790  7.674084e-05 1.0003675
## prev[2]   1.249699e-04     1.0 10836 -3.169610e-04 1.0003569
## prob_1[1] 6.537476e-05     1.0 10905 -4.468137e-04 1.0003389
## prob_1[2] 6.259661e-05     1.1  8063 -2.809511e-04 1.0001406
## prob_1[3] 4.322508e-05     1.1  8433  7.642535e-03 1.0002174
## prob_1[4] 9.264131e-05     1.0 10541  5.740739e-04 1.0000597
## prob_2[1] 1.187534e-04     0.9 12095  4.492070e-03 1.0003792
## prob_2[2] 4.101602e-05     0.9 12366  1.028519e-03 0.9999808
## prob_2[3] 4.206362e-05     1.0 10418  2.047872e-02 1.0000397
## prob_2[4] 1.149664e-04     0.9 11707  1.745897e-03 1.0002765
## se[1]     6.921380e-05     1.1  8727  2.335653e-02 1.0003048
## se[2]     7.266383e-05     1.0  9427  7.800848e-03 0.9999859
## sp[1]     6.640189e-05     1.1  7935  1.804983e-04 1.0001219
## sp[2]     4.661818e-05     1.1  8209  7.924786e-03 1.0003453

Exercise 3

Try to run the model with different initial values, that explore the whole parameter space and removing I(1-sp[1],).

For example try it with:

se <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))
sp <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))

Solution

The model is the same like as in the above scenario (hw_definition_2)

Population_1 <- as.numeric(pop_1[1:2,1:2])
Population_2 <- as.numeric(pop_2[1:2,1:2])


prev <- list(chain1=c(0.05,0.99), chain2=c(0.95,0.05))
se <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))
sp <- list(chain1=c(0.01,0.99), chain2=c(0.99,0.01))

results_3 <- run.jags(hw_definition_2, n.chains=2, sample = 100000)
# Remember to check convergence and effective sample size!
plot(results_3)  # Click backwards to view all plots

pt_3 <- plot(results_3)

pt_3$`prev[1].plot1`

pt_3$`prev[1].plot3`

print(pt_3[["prev[1].plot1"]])

print(pt_3[["se[1].plot1"]])

print(pt_3[["sp[1].plot1"]])

print(pt_3[["sp[1].plot3"]])

print(pt_3[["se[2].plot1"]])

summary(results_3)
##               Lower95      Median    Upper95        Mean          SD Mode
## prev[1]   0.015154582 0.028079649 0.04294688 0.028632057 0.007270327   NA
## prev[2]   0.690483159 0.715630958 0.74086809 0.715340989 0.012949846   NA
## prob_1[1] 0.014662730 0.026397736 0.04058197 0.026919481 0.006779520   NA
## prob_1[2] 0.007950100 0.017411169 0.02927751 0.018002716 0.005609978   NA
## prob_1[3] 0.002160442 0.008386046 0.01666301 0.008947818 0.003947023   NA
## prob_1[4] 0.927542366 0.946665368 0.96424735 0.946129985 0.009514629   NA
## prob_2[1] 0.642836937 0.669039765 0.69347360 0.669034858 0.012993111   NA
## prob_2[2] 0.020017827 0.028496244 0.03771279 0.028727058 0.004575018   NA
## prob_2[3] 0.015947471 0.023963876 0.03237984 0.024224011 0.004216653   NA
## prob_2[4] 0.253945453 0.277856590 0.30219512 0.278014073 0.012380665   NA
## se[1]     0.956013612 0.968677497 0.98071913 0.968372692 0.006350692   NA
## se[2]     0.951911789 0.965938232 0.97953781 0.965773842 0.007061216   NA
## sp[1]     0.970692441 0.982929910 0.99304449 0.982298017 0.005902372   NA
## sp[2]     0.983253015 0.992136762 0.99871786 0.991542641 0.004195551   NA
##                  MCerr MC%ofSD SSeff         AC.10      psrf
## prev[1]   5.140897e-05     0.7 20000  2.851607e-04 1.0004309
## prev[2]   9.222759e-05     0.7 19715  3.357144e-03 1.0000006
## prob_1[1] 4.793845e-05     0.7 20000  2.100863e-04 1.0004153
## prob_1[2] 4.006850e-05     0.7 19603  1.397071e-03 0.9999973
## prob_1[3] 2.814228e-05     0.7 19671 -1.038447e-03 1.0001594
## prob_1[4] 6.727859e-05     0.7 20000  2.223914e-04 1.0002099
## prob_2[1] 9.187517e-05     0.7 20000 -5.566153e-04 1.0000736
## prob_2[2] 3.235026e-05     0.7 20000  8.991125e-04 0.9999931
## prob_2[3] 2.981624e-05     0.7 20000 -8.326157e-04 1.0003006
## prob_2[4] 8.754452e-05     0.7 20000  2.277412e-03 1.0000170
## se[1]     4.490618e-05     0.7 20000 -3.803709e-05 1.0002579
## se[2]     4.993034e-05     0.7 20000  4.214278e-03 0.9999810
## sp[1]     4.214329e-05     0.7 19615  1.900281e-03 1.0000090
## sp[2]     2.992987e-05     0.7 19650 -8.217058e-04 1.0003295

Exercise 4

Run the model with only 1 population (either pop_1 or pop_2). What happens then?

Solution

Model Specification (‘hw_definition’)

hw_definition_4 <- c("model{
  Population_1 ~ dmulti(prob_1, N_1)

  #Population_1
  
  # Test1+ Test2+
    prob_1[1] <- (prev[1] * ((se[1])*(se[2]))) + ((1-prev[1]) * ((1-sp[1])*(1-sp[2])))
  
  # Test1+ Test2-
    prob_1[2] <- (prev[1] * ((se[1])*(1-se[2]))) + ((1-prev[1]) * ((1-sp[1])*(sp[2])))

  # Test1- Test2+
    prob_1[3] <- (prev[1] * ((1-se[1])*(se[2]))) + ((1-prev[1]) * ((sp[1])*(1-sp[2])))

  # Test1- Test2-
    prob_1[4] <- (prev[1] * ((1-se[1])*(1-se[2]))) + ((1-prev[1]) * ((sp[1])*(sp[2])))
    

  prev[1] ~ dbeta(1, 1)
  se[1] ~ dbeta(1, 1)I(1-sp[1],)
  sp[1] ~ dbeta(1, 1)
  se[2] ~ dbeta(1, 1)I(1-sp[2],)
  sp[2] ~ dbeta(1, 1)

  #data# Population_1, N_1
  #monitor# prev, prob_1,  se, sp
  #inits# prev, se, sp
  }
  ")

Population_1 <- as.numeric(pop_1[1:2,1:2])


prev <- list(chain1=c(0.05), chain2=c(0.95))
se <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))
sp <- list(chain1=c(0.5,0.99), chain2=c(0.99,0.5))

results_4 <- run.jags(hw_definition_4, n.chains=2)
# Remember to check convergence and effective sample size!
plot(results_4)  # Click backwards to view all plots

pt_4 <- plot(results_4)

pt_4$`prev.plot1`

pt_4$`prev.plot3`

print(pt_4[["prev.plot1"]])

print(pt_4[["se[1].plot1"]])

print(pt_4[["sp[1].plot1"]])

print(pt_4[["sp[1].plot3"]])

print(pt_4[["se[2].plot3"]])

print(pt_4[["se[2].plot4"]])

print(pt_4[["se[2].plot5"]])

summary(results_4)
##               Lower95      Median    Upper95       Mean          SD Mode
## prev      0.019013278 0.038341165 0.06161965 0.03946337 0.011228794   NA
## prob_1[1] 0.012043887 0.023475571 0.03662748 0.02406689 0.006416080   NA
## prob_1[2] 0.008585623 0.018714963 0.03110883 0.01933945 0.005895514   NA
## prob_1[3] 0.003211859 0.009860699 0.01886579 0.01044148 0.004226694   NA
## prob_1[4] 0.927745026 0.946748161 0.96475071 0.94615218 0.009599461   NA
## se[1]     0.634003775 0.852750996 0.99999306 0.83821366 0.109583867   NA
## se[2]     0.506944192 0.754465929 0.99999240 0.75081035 0.144351054   NA
## sp[1]     0.976337501 0.989594679 0.99999630 0.98892152 0.007029213   NA
## sp[2]     0.986310908 0.994833454 0.99999952 0.99411332 0.004161227   NA
##                  MCerr MC%ofSD SSeff        AC.10      psrf
## prev      1.661421e-04     1.5  4568  0.022368201 1.0003222
## prob_1[1] 5.537503e-05     0.9 13425 -0.013294259 1.0004676
## prob_1[2] 5.279877e-05     0.9 12468  0.002956750 1.0000889
## prob_1[3] 4.239210e-05     1.0  9941  0.002543146 1.0001484
## prob_1[4] 8.900707e-05     0.9 11632 -0.001391273 1.0000299
## se[1]     1.471751e-03     1.3  5544 -0.008164718 0.9999978
## se[2]     2.344950e-03     1.6  3789  0.045963785 1.0002313
## sp[1]     1.097099e-04     1.6  4105  0.041163880 1.0001345
## sp[2]     5.603503e-05     1.3  5515 -0.002938294 1.0000742

How do the results look? # Compare the results with the scenarios above? Try also by removing population 1